function [p,e,t] = grid_ff2oct(gridfile)
% [p,e,t] = grid_ff2oct(gridfile)
%
% This function reads a FreeFem++ grid in format .mesh from file
% "gridfile" and reformats it in the usual [p,e,t] "Matlab-like"
% format.
% LIMITATIONS: 
% - only one label for the domains is recognized
% - the curvilinear abscissa in e is not set

 fid = fopen(gridfile,'r');

 % determine the dimension
 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   pos = strfind(line,'Dimension');
 end
 if(length(line)>length('Dimension')+1) % line includes a \n
   line(pos+length('Dimension'):end);
   dim = sscanf(line(pos+length('Dimension'):end),'%d',1);
 else
   line = fgets(fid);
   dim = sscanf(line,'%d',1);
 end

 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   pos = strfind(line,'Vertices');
 end
 nv = fscanf(fid,'%d',1);
 % vertex coordinates
 if(dim==2)
   xyz = fscanf(fid,'%f %f %d',[3,nv]);
 else
   xyz = fscanf(fid,'%f %f %f %d',[4,nv]);
 end
 % discard the label
 p = xyz(1:dim,:);

 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   if(dim==2)
     pos = strfind(line,'Edges');
   else
     pos = strfind(line,'Triangles');
   end
 end
 ns = fscanf(fid,'%d',1);
 % sides
 if(dim==2)
   xyz = fscanf(fid,'%d %d %d',[3,ns]);
 else
   xyz = fscanf(fid,'%d %d %d %d',[4,ns]);
 end
 e = zeros(dim+3,ns);
 e(1:dim,:) = xyz(1:dim,:);
 e(dim+3,:) = xyz(dim+1,:);

 fseek(fid,0);
 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   if(dim==2)
     pos = strfind(line,'Triangles');
   else
     pos = strfind(line,'Tetrahedra');
   end
 end
 ne = fscanf(fid,'%d',1);
 % elements
 if(dim==2)
   xyz = fscanf(fid,'%d %d %d %d',[4,ne]);
 else
   xyz = fscanf(fid,'%d %d %d %d %d',[5,ne]);
 end
 % discard the label
 t = xyz(1:dim+1,:);

 fclose(fid);

return

